
Generate_Contract_Specs <- function(input,
                              r = 0.03,
                              lambda = 0.20,
                              mu_d_RN = 0.0735,
                              sigma_d = 0.14,
                              mu_j = -0.25,
                              sigma_j = 0.10,
                              Cost_grid = c(0.002, 0.003, 0.005, 0.007),
                              N_sim = 1e+6,
                              MRP = 0.06)
{
  
  # Testing
  # input <- list()
  # input$age = 35
  # input$retirement_age = 65
  # input$income = 1000
  # 
  # 
  # r = 0.03 # constant annual interest rate, compounded continuously.
  # lambda = 0.20 # likelihood of jump.
  # mu_d_RN = 0.0735 # mean of diffusion process (under risk-neutral measure).
  # sigma_d = 0.14 # volatility of diffusion process.
  # mu_j = -0.25 # mean of each jump.
  # sigma_j = 0.10 # volatility of magnitude of each jump.
  # Cost_grid = c(0.002, 0.003, 0.005, 0.007) # cost ('fee_TDF' or 'c_bar_RILA').
  # N_sim = 1e+4 # number of simulated paths.
  # MRP = 0.06 # Market risk premium
  
  
  
  
  
  
  # Vars from input
  age <- input$age
  retirement_age <- input$retirement_age
  income <- input$income  
  
  # Other vars
  mu_d = mu_d_RN + MRP # mean of diffusion process (under real-world measure).
  C_num = length(Cost_grid)  
  
  
  # Function to help find fair RILA cap rates -------------------------------
  ExpLoss <- function(Cap,
                      D_level,
                      D_type,
                      r,
                      tau,
                      fee,
                      c_bar,
                      R_S,
                      A_t = 1000)
  {
    # Cap = 1
    # D_level = con_spec_vec[t,j]
    # D_type = "Buffer"
    # r = r
    # tau = 1
    # fee = 0
    # c_bar = c_bar_RILA
    # logR = logR
    A_t = 1000
    
    # simulated annual net returns of market index.
    A_t = 1000
    
    if(D_type == "Buffer"){
      
      # vector (N_sim-by-1) of credited return to RILA account.
      R_A = pmin(R_S + D_level, 0) * (R_S < 0) + pmin(R_S, Cap) * (R_S > 0) 
      
    } else if(D_type == "Floor"){
      
      # vector (N_sim-by-1) of credited return to RILA account.
      R_A = pmax(R_S , -D_level) * (R_S < 0) + pmin(R_S, Cap) * (R_S > 0) 
      
    } else {
      print('Error. Downside protection type unspecified.')
      asdfadsf    # to produce error
    }
    
    A_tplustau = A_t * exp(-fee * tau) * (1 + R_A)
    L = A_tplustau * exp(-r * tau) + A_t * c_bar * tau - A_t
    EL = mean(L)
    
    return(EL)
  }
  
  # Contract specs ----------------------------------------------------------
  
  # number of contracts used in experiment.
  num_con = 9                    
  
  con_type_vec <- c()
  con_type_vec[1:3] = "TDF"
  con_type_vec[4:6] = "Buffer"
  con_type_vec[7:9] = "Floor"
  
  # contains 'alpha' for TDF, 'B' for buffer-RILA and 'F' for floor-RILA.
  con_spec_vec = matrix(0, 50, num_con)
  
  # TDFs:
  con_spec_vec[1:20, 1] = 1       
  con_spec_vec[21:50, 1] = 0.95 - 0.01 * c(1:30)
  con_spec_vec[, 2] = con_spec_vec[, 1] - 0.10
  con_spec_vec[, 3] = con_spec_vec[, 1] - 0.20
  
  # Buffer RILAs:
  con_spec_vec[1:20, 4] = 0     
  con_spec_vec[21:50, 4] = 1/300 * c(1:30)
  con_spec_vec[, 5] = con_spec_vec[, 4] + 0.025
  con_spec_vec[, 6] = con_spec_vec[, 4] + 0.050
  
  # Floor RILAs:
  con_spec_vec[1:20, 7] = 0.25  	
  con_spec_vec[21:50, 7] = 0.25 - 1.9/300 * c(1:30)
  con_spec_vec[, 8] = con_spec_vec[, 7] - 0.035
  con_spec_vec[, 9] = con_spec_vec[, 7] - 0.055
  
  
  # Find Cap Rates of RILA contracts: ---------------------------------------
  
  # contains cap rates for RILA contracts.
  con_cap_vec = array(0, dim = c(50, C_num, num_con))     
  
  # simulate log returns of index for 1 year:
  num_jumps = rpois(N_sim, lambda)
  
  # vector of annual log returns, following MJD model.
  logR = (mu_d - sigma_d^2/2)*1 + 
    num_jumps * mu_j + 
    sqrt( sigma_d^2/2 * 1 + num_jumps * sigma_j^2 ) * rnorm(n = N_sim, mean = 0, sd = 1)   
  
  R_S = exp(logR) - 1 
  
  
  con_spec_vec_comb <- con_spec_vec[t,j]
  for(t in 1:50){
    
    for(c in 1:C_num){
      
      for(j in 1:num_con){
        
        if(con_type_vec[j] == "TDF"){
          # do nothing
          
        } else if(con_type_vec[j] == "Buffer") {
          con_cap_vec[t,c,j] = uniroot(ExpLoss, 
                                       interval = c(0, 1),
                                       D_level = con_spec_vec[t,j],
                                       D_type = "Buffer",
                                       r = r,
                                       tau = 1,
                                       fee = 0, 
                                       c_bar = Cost_grid[c],
                                       R_S = R_S,
                                       tol = .01)$root
          
          
          
        } else if(con_type_vec[j] == "Floor") {
          con_cap_vec[t,c,j] = uniroot(ExpLoss, 
                                       interval = c(0, 1),
                                       D_level = con_spec_vec[t,j],
                                       D_type = "Floor",
                                       r = r,
                                       tau = 1,
                                       fee = 0, 
                                       c_bar = Cost_grid[c],
                                       R_S = R_S,
                                       tol = .01)$root
          
        } else {
          print('Contract type unspecified.')
          
        }
        
      }
      print(t*c*j/(50 * C_num * num_con))
    }
    
  }
  
  output <- list(
    Cost_grid = Cost_grid,
    con_type_vec = con_type_vec,
    con_spec_vec = con_spec_vec,
    con_cap_vec = con_cap_vec)
  
  return(output)
}